NASA Computational Case Study 
The Flight of Friendship 7 


David G. Simpson 

NASA Goddard Space Flight Center, Greenbelt, Maryland, 20771 


NASA Science Application: Astrodynamics, celestial mechanics, naviga- 
tion, Project Mercury 

Computational Algorithms: Modeling, time and coordinate conversions, 
root finding 


Abstract 

In this case study, we learn how to compute the position of an Earth-orbiting 
spacecraft as a function of time. As an exercise, we compute the position of 
John Glenn’s Mercury spacecraft Friendship 7 as it orbited the Earth during 
the third flight of NASA’s Mercury program. 

1 Introduction 

On February 20, 1962, astronaut John Glenn became the third American as- 
tronaut in space and the first American to orbit the Earth in his spacecraft 
Friendship 7. Glenn’s spacecraft was launched into orbit atop a Mercury- Atlas 
rocket at 9:47:39 am EST from Cape Canaveral, Florida (Figure 1). He com- 
pleted three orbits around the Earth in just under five hours, after which his 
spacecraft splashed down in the Atlantic Ocean at 2:43:09 pm EST, about 260 
miles north of San Juan, Puerto Rico. He was picked up by the destroyer USS 
Noa a few minutes after splashdown. 

Friendship 7, like any Earth-orbiting spacecraft, traveled in an elliptical orbit 
about the Earth, with the Earth’s center of mass at one focus of the ellipse. The 
point in the orbit closest to Earth is called perigee ; the point farthest from Earth 
is called apogee. For Friendship 7, the perigee altitude was 86.92 nautical miles 1 
above the Earth’s surface, and the apogee altitude was 140.92 nautical miles. [1] 
Calculating the position of Friendship 7 at any time involves the application 
of some basic methods of celestial mechanics, which we’ll examine here. 

X A un ' t °f distance frequently used in spacecraft work: 1 nautical mile is exactly 1852 
meters. 
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Figure 1: Launch of the Mercury-Atlas 6 mission on February 20, 1962. The 
Friendship 7 spacecraft is the dark cone at the top of the Mercury- Atlas rocket 
(Credit: NASA ) 

2 Time Measurement 

In order to calculate the position of an Earth-orbiting spacecraft, we begin 
with finding a suitable method for measuring time. In everyday civil use, we 
measure time using a rather cumbersome system: years, months (of irregular 
length, between 28 and 31 days), and day of month. Within a single day, we 
divide time into hours, minutes, and seconds. 

This system is not particularly convenient for use in calculations or plotting, 
where we would like to have time varying smoothly along some continuous time 
scale and measured with a single unit. Astronomers have developed such a 
uniform scale of time measurement, called the Julian day: it is defined to be the 
total number of days elapsed since noon on Monday, January 1, 4713 B.C. [2,3] 
For example, noon on January 1, 2000 is Julian day 2451545.0, since that’s the 
number of days that had elapsed in the 6712 years since the beginning of 4713 
B.C. 

A fractional day may be added to the Julian day. Usually the fractional 
day is counted from midnight Coordinated Universal Time (UTC), which is five 
hours ahead of Eastern Standard Time (EST). 2 

The Julian day for any date on the Gregorian calendar may be found by 
consulting a table of Julian days [4], or computed using a simple algorithm. Let 
Y by the year M the month number (1 for January, 2 for February, etc., up to 12 
for December), and D be the day of month (including a fractional day, counted 
from midnight). The algorithm for computing the Julian day (JD) is [3]: 

2 Other time scales may be used as well. When it’s important to be clear which time scale 
is used in computing the fractional day, the time scale is added in parentheses following JD: 
for example, “JD(UTC)”. 
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• If M > 2, then leave Y and M unchanged. 

If M = 1 or M = 2, then replace Y with Y — 1, and M with M + 12. 

• Calculate 

^ = INT (lo) (') 

B = 2-A + INt(+) (2) 

• Then the Julian day JD is found from 

JD = INT[365.25(T +4716)] +INT[30.6001(Af +1)] +D + B- 1524.5 (3) 
Here INT(x) indicates the greatest integer less than or equal to x. 


EXAMPLE 1. Robert Goddard launched his first experimental rocket in Auburn, 
Massachusetts, on March 16, 1926, at 19:30 UTC. What w^as the corresponding 
Julian day? 


Solution. The time 19:30 corresponds to a fractional day 
19 30 

o7 TTTrv ~ 0.8125 day, 

24 1440 

where we have used 1 day = 24 hours = 1440 minutes. Now using the above 
algorithm to compute the Julian day, we find: 

Y = 1926, M = 3, D = 16.8125, A = 19, B = —13, 

and so 


JD - 2425990 + 122 + 16.8125 


13 — 1524.5 = 


2424591.3125 JD 


Activity 1 . Apollo 11 astronaut Neil Armstrong first set foot on the Moon on 
July 21, 1969, at 02:56 UTC. Find the Julian day corresponding to this time. 


3 Reference Frames 

In order to describe an orbit mathematically, it is necessary to introduce a refer- 
ence frame, with respect to which the orientation of the orbit can be measured. 
Such a reference frame is determined by a reference plane, as well as a reference 
direction within that plane. Two reference planes are in common use: 
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• The equator is the plane of the Earth’s equator. This is the reference 
plane usually chosen for bodies orbiting the Earth. 

• The ecliptic is the plane of the Earth’s orbit. This is the reference plane 
chosen in most other cases: finding positions of Sun-orbiting bodies such 
as planets, comets, asteroids, etc. 3 

In both cases, the reference direction is chosen to be the direction of the vernal 
equinox , which is in the direction from the Earth to the Sun on the first day of 
spring. 4 It is a fixed point in space in the constellation Pisces, and lies along 
the line of intersection between the planes of the equator and the ecliptic. The 
line pointing toward the vernal equinox is therefore common to both reference 
planes. 

When an orbiting body crosses the reference plane from south to north , it 
is said to be at the ascending node of the orbit. Crossing the plane in the other 
direction, north to south, is the descending node. 

Throughout this case study, we’ll use the equator as the reference plane. 


4 Orbital Elements 

A spacecraft orbiting the Earth will travel in an elliptical orbit, with the center 
of mass of the Earth at one focus of the ellipse. To calculate the position of an 
orbiting body at any time, we need several parameters to describe the orbit: 

• The size of the orbit is given by the semi-major axis 5 a of the ellipse. 

• The shape of the orbit is given by the eccentricity e of the ellipse. For an 
ellipse with semi-major axis a and semi-minor axis 6, the eccentricity is 
d efined to be the ratio of the distance between the center and either focus 
\/ a 2 — b 2 to the semi-major axis a: 


\/a 2 — b 2 

e = 

a 


(4) 


The eccentricity is a dimensionless number that varies between 0 (for a 
perfect circle) to almost 1 (for a long, elongated, cigar-shaped ellipse). 

• The orientation of the orbit in space with respect to the equatorial refer- 
ence frame is determined by three angles, as shown in Figure -2: 

- The inclination i of the orbit. This is the dihedral angle between the 
orbit plane and the plane of the Earth’s equator. 

3 Although it orbits the Earth, the position of the Moon is usually described using the 
ecliptic as the reference plane. 

4 That is, spring in the northern hemisphere — around March 21. 

5 The semi-major axis a is half the long axis of the ellipse; the semi-minor axis b is half 
the short axis. 
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Figure 2: Orbital elements. In the case of an Earth-orbiting body, the xy plane 
is the plane of the equator, x is in the direction of the vernal equinox, and the 
2 axis points northward along the Earth’s rotation axis. Here the Earth is at 
the origin, and m is the orbiting body; i is the inclination, Q is the longitude of 
the ascending node IV, lu is the argument of perigee, and / is the true anomaly. 
(After [5].) 


— The longitude of the ascending node fl. This is the angle, measured 
in the plane of the Earth’s equator, between the vernal equinox and 
the ascending node of the orbit. 

- The argument of perigee to. This is the angle, measured in the plane 
of the orbit, between the ascending node and the perigee point. 

• The five elements given so far completely describe the orbit in space. We 
now only need to specify where in the orbit the spacecraft can be found 
at some given epoch time T 0 . This is given by an angle called the mean 
anomaly at the epoch time , M 0 . Angle M 0 is the angle, measured in the 
plane of the orbit, from the perigee point to the spacecraft position at 
epoch time T 0 , assuming the spacecraft moves at a constant rate around 
the Earth. 


5 Computation of the Position 

We can now begin the orbit calculations, for which we’ll use standard two- 
body orbit propagation methods. [5] Our goal will be to find the latitude and 
longitude L of an Earth-orbiting spacecraft for a given time t. 
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First, we will need to find the mean daily motion n of the orbit, which is 
just the number of revolutions (orbits) the spacecraft makes around the Earth 
per day. The mean daily motion may be given, or it may be computed from 
Kepler’s third law if the semi-major axis a is known: 


_ 86400 j GM e 

2ir V a 3 


( 5 ) 


Here G is Newton’s gravitational constant, M e is the mass of the Earth, a is the 
semi-major axis of the orbit, and the factor 86400/27T converts from SI units of 
rad/s to units of rev/day. The product GM e has a value 6 ’ 7 of 3.986004415 x 10 14 
m 3 s -2 . 

Next, we begin by assuming the spacecraft is moving in a circular orbit at 
a constant rate; we’ll correct for the eccentricity of the orbit later. Assuming 
a circular orbit, the mean anomaly M of the spacecraft at time t can be found 
from knowing the mean anomaly M 0 at the epoch time T 0 : 


Af — Mo - J - 2tt n(t — To). (6) 

Here M is the mean anomaly at time t, M 0 is the mean anomaly at the epoch 
time T 0 , and n is the mean daily motion. Both M and M 0 are in units of 
radians , both t and To are specified as Julian days, and n has units of rev/day. 
(The factor of 2n converts the second term on the right from units of revolutions 
to radians.) 

Next we introduce a correction for the eccentricity of the orbit. We do this 
by computing an angle called the eccentric anomaly E, which is found by solving 
Kepler’s equation , 

M — E — esinT. (7) 

Here both M and E must be in radians, 8 and e is the given eccentricity. The 
mean anomaly M is known from Eq. (6), and we wish to solve for the eccentric 
anomaly E. However, Kepler’s equation cannot be solved explicitly for E\ we 
must instead employ a numerical method to solve it. This is discussed in detail 
later. 

The next step is to find another angle, the true anomaly /. This is the angle, 
measured in the plane of the orbit, from the perigee point to the spacecraft’s 
position at time t, with the focus of the orbit at the vertex of the angle (Figure 
2). The true anomaly / is found from 

(£) = (rrf ) 7 tan (!) ■ (*) 

6 As given in Ref. [4]. The product GM e is known to better accuracy than either G or M e 
individually, so it’s best to use this product in computing Eq. (5). 

7 The reciprocal of n is the orbital period , or time required to orbit the Earth, in days. 

When an equation contains a “bare” angle that’s not the argument of a trigonometric 
function as M and the first E are here — the angle is assumed to be in radians. 
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We now need to compute the radial distance r from the orbit focus to the 
spacecraft. This is found from the eccentric anomaly E using 

r = a( 1 — ecos E), (9) 

where r will have the same units as the given semi-major axis a. 

The quantities (r, /) are the polar coordinates of the spacecraft, in the plane 
of the orbit; the rest of the work is a series coordinate transformations, to go from 
plane polar coordinates in the plane of the orbit to spherical polar coordinates 
with the xy plane at the equator. We begin these coordinate transformations 
by defining the argument of latitude u by 

u — ^ + /? ( 10 ) 

where is the given argument of perigee and / is the true anomaly. Then 
the cartesian coordinates of the spacecraft in the equatorial frame at time t are 
given by 


X 

= r (cos u cos Vt — sin u sin Q cos i) 

(ii) 

y 

= r (cos u sin O + sin u cos Q cos i) 

(12) 

z 

= r sin u sin i 

(13) 


Converting from cartesian to spherical polar coordinates gives an azimuthal 
angle a called the right ascension , and a co-polar angle 5 called the declination 
(Figure 3): 


V 

tan a = — 

x 

sind = - 

r 


(14) 

(15) 


Finally, the longitude L and latitude ip are related to the right ascension and 
declination through a rotation about the z axis due to the Earth’s rotation: 


L = a- GST (16) 

<p = 5 (17) 

Here GST is the Greenwich Sidereal Time, and is the angle from the vernal 
equinox to the prime meridian (measured in the plane of the equator) at time 
t. The GST (in degrees) may be found from an empirical formula [3]: 

GST = 280.46061837 + 360.98564736629(t — 2451545.0) 

+ 0.000387933 T 2 - T 3 /38710000, (18) 

where t is a Julian day, and T is the time in Julian centuries (of 36525 days) 
from noon, January 1, 2000: 


t - 2451545.0 
36525 


( 19 ) 
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Figure 3: Right ascension a and declination 5 of an object, shown here as a 
star. As in Fig. 2, the origin is at the center of the Earth, the x axis is toward 
the vernal equinox, and the z axis points northward along the Earth’s rotation 
axis. (After [6].) 


Equation (18) may return very large angles; you should add or subtract multiples 
of 360° as needed to bring the GST into the range [ 0°, 360°). 

Once you’ve computed the longitude L using Eq. (16), you should reduce it to 
the range [-180°, +180°] by adding or subtracting multiples of 360° as needed. 
Then positive L corresponds to east longitude, and negative L to west longitude. 
The latitude p should already be in the range [-90°, +90°] after taking the 
inverse sine of both sides of Eq. (15); in this case positive ip corresponds to 
north latitude, and negative ip to south latitude. 


Activity 2. The first satellite launched by the United States was Explorer 1, 
a 30-lb satellite in the shape of a cylinder about 5 feet long and 6 inches in 
diameter. It was launched into orbit by a Jupiter-C rocket on January 31, 
1958, and was responsible for the discovery of the Van Allen radiation belts. It 
remained in orbit until 1970. [7] 

The orbit of Explorer 1 had a semi-major axis a of 7615.480 km, and an 
eccentricity e of 0.1155556. [8] (a) What was the mean daily motion n? ( b ) What 
was the orbital period? (c) What were the altitudes of apogee and perigee? 
(Assume for simplicity that the Earth is a sphere of radius 6378 km.) 
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POINTER. Inverse Trigonometric Functions. 

Exercise care when computing inverse trigonometric functions, to ensure that 
the resulting angle is in the correct quadrant. In general, there are two correct 
angles, but the inverse trigonometric function of a calculator or computer will 
return only one of them. The calculator-provided angle will be between —90° and 
+90° for the sin -1 and tan -1 functions, and between 0° and 180° for the cos -1 
function. 

In a case like Eq. (14), where we compute the inverse tangent of the ratio 
y/x, determining the correct quadrant is simple: if the denominator x is negative , 
then add 180° to the calculator’s or computer’s returned answer; otherwise use 
the answer as returned. Many computer programming languages provide a special 
function (typically called something like atan2 (y, x ) ) just for this type of problem, 
which will automatically return the angle in the correct quadrant. 

In the case of Eq. (15), the returned angle will be between -90° and +90°; 
but since 5 is a “latitude” angle, that’s where it should be, so no adjustment of 
the angle is necessary. 

Eq. (8) is a bit more complicated to analyze, but it turns out that / will 
automatically be in the correct quadrant using the returned values for the tan -1 
function. 

Remember that you can always add or subtract as many multiples of 360° as 
i you like without changing the angle. 


Activity 3. (a) Referring to your solution to Activity 1, compute the Green- 
wich Sidereal Time (GST) for the moment Neil Armstrong first stepped onto 
the Moon. Reduce your answer to the range [0°,360°). (b) Given that the 
right ascension of the Moon at the time was 190.2°, what was the geographic 
longitude of the Moon? (c) What part of the Earth would have been visible to 
Armstrong when he first looked back at the Earth from the Moon? 


6 Solving Kepler’s Equation 

The calculations just described involved solving Kepler’s equation, which relates 
the mean anomaly M to the eccentric anomaly E : 


M = E — e sin E, (20) 

where e is the eccentricity of the orbit, and both M and E must be in radians. 
We are given M and e, and wish to solve for E. This cannot be done in closed 
form, but must be done numerically. 

Quite a few numerical methods for solving Kepler’s equation have been de- 
veloped over the years. [9] One simple method of solution is Newton’s method. 
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Given a function /(x), Newton’s method is an iterative technique for finding 
the roots of /(x), i.e. those values of x for which 

/(*) = 0 . ( 21 ) 

To find a root using Newton’s method, we begin with an initial estimate xo of 
the root. Then successively better approximations to the root are found using 
the iterative formula [10] 

f(Xi) 

X i+ 1 ~ x i 777 7 5 (22) 

where f'(x) is the first derivative of /(x). Starting with the initial estimate 
xo, this formula gives a better estimate x\. Substituting this x\ back into 
the formula again gives an even better estimate x 2 , and so on. Repeating this 
process gives better and better estimates of the root of /(x); you simply continue 
the process until the answer converges to the desired accuracy (i.e. the difference 
between the current and previous iterations is below a certain threshold). 

In the case of Kepler’s equation, we wish to find the root E of the equation 

f(E) = M — E + e sin E = 0. (23) 

Therefore the iteration formula for the solution to Kepler’s equation using New- 
ton’s method becomes 


E, 


n+l 


En. ~ 


M — E n + e sin E n 
e cos E n — 1 


(24) 


Newton’s method requires an initial estimate for E, for which you can just use 
the mean anomaly M, so that Eq = M. 


Activity 4. Derive Eq. (24) from Eqs. (22) and (23). 


Activity 5. Halley’s comet is in a highly elongated elliptical orbit, with eccen- 
tricity e = 0.967. If the mean anomaly of Halley’s comet is M = 215°, then use 
Newton’s method to solve Kepler’s equation for the eccentric anomaly E. 

Hint 1: You will need to work this problem in radians. Begin by converting M 
from degrees to radians (by multiplying by 7r/180). 

Hint 2: Use M as the initial estimate of E, then apply Eq. (24) repeatedly. It 
should only take 3 iterations for the answer to converge to 1 1 significant digits. 
Hint 3: At the end, convert your final answer E from radians back to degrees 
(by multiplying by 180 /tt). 


10 


Challenge. Now let’s compute the position of the Friendship 7 spacecraft at a 
specific time. During his first orbit, John Glenn reported seeing what he called 
“fireflies” — small lighted objects swirling outside the spacecraft: 

“I am in a big mass of some very small particles, that are brilliantly 
lit up like they’re luminescent. I never saw anything like it. They 
round a little; they’re coming by the capsule, and they look like 
little stars. A whole shower of them coming by. They swirl around 
the capsule and go in front of the window and they’re all brilliantly 
lighted. They probably average maybe 7 or 8 feet apart, but I can 
see them all down below me, also.” 

It was later determined that the most likely cause of these lights was bits of 
paint or other material flaking off the spacecraft. 

Glenn reported seeing these “fireflies” at around time 16:03:03 UTC. Find 
the latitude and longitude of the Friendship 7 spacecraft at this time, given the 
orbital elements of the spacecraft in Table 1 below. 


Table 1. Orbital elements of Friendship 7. 


Orbital element 

Symbol 

Value 

Semi-major axis 

a 

6589.116 km 

Eccentricity 

e 

0.007589 

Inclination 

i 

32.54° 

Longitude of ascending node 

0 

235.2° 

Argument of perigee 

L J 

181.2° 

Mean anomaly at epoch 

Mo 

228.5° 

Epoch time 

To 

2437716.11642 JD 


If you write a computer program to compute the positions of the spacecraft 
at, for example, 10-second intervals during the entire mission — from launch to 
splashdown — and plot them on a map of the Earth, you’ll see the sinusoidal 
ground track of Friendship 7, which shifts westward each orbit due to the rota- 
tion of the Earth. 
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